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Context. The Tarantula Nebula in the Large Magellanic Cloud is our closest view of a starburst region and is the ideal environment to investigate 
important questions regarding the formation, evolution and final fate of the most massive stars. 
Q_l Aims. We analyze the multiplicity properties of the massive O-type star population observed through multi-epoch spectroscopy in the framework 
D of the VLT-FLAMES Tarantula Survey. With 360 O-type stars, this is the largest homogeneous sample of massive stars analyzed to date. 

Methods. We use multi-epoch spectroscopy and variability analysis to identify spectroscopic binaries. We also use a Monte-Carlo method to 
correct for observational biases. By modeling simultaneously the observed binary fraction, the distributions of the amplitudes of the radial velocity 
variations and the distribution of the time scales of these variations, we derive the intrinsic current binary fraction and period and mass-ratio 
distributions. 

Results. We observe a spectroscopic binary fraction of 0.35 ± 0.03, which corresponds to the fraction of objects displaying statistically significant 
radial velocity variations with an amplitude of at least 20 km s -1 . We compute the intrinsic binary fraction to be 0.5 1 ± 0.04. We adopt power-laws 
to describe the intrinsic period and mass-ratio distributions: /(log 10 P/d) ~ (log 10 P/d)" (with log 10 P/d in the range 0.15-3.5) and f(q) ~ q" with 
0.1 < q = M 2 /M[ < 1.0. The power-law indexes that best reproduce the observed quantities are n = -0.45 ± 0.30 and k = -1.0 ± 0.4. The period 
r~ | distribution that we obtain thus favours shorter period systems compared to an Opik law (n = 0). The mass ratio distribution is slightly skewed 
O | towards low mass ratio systems but remains incompatible with a random sampling of a classical mass function (k = -2.35). The binary fraction 
I seems mostly uniform across the field of view and independent of the spectral types and luminosity classes. The binary fraction in the outer 
O region of the field of view (r > 7.8', i.e. «1 17 pc) and among the 09.7 I/II objects are however significantly lower than expected from statistical 
fluctuations. The observed and intrinsic binary fractions are also lower for the faintest objects in our sample (K s > 15.5 mag), which results from 
C/J observational effects and the fact that our O star sample is not magnitude-limited but is defined by a spectral-type cutoff. We also conclude that 
magnitude-limited investigations are biased towards larger binary fractions. 

Conclusions. Using the multiplicity properties of the O stars in the Tarantula region and simple evolutionary considerations, we estimate that over 
, 50% of the current O star population will exchange mass with its companion within a binary system. This shows that binary interaction is greatly 
^ affecting the evolution and fate of massive stars, and must be taken into account to correctly interpret unresolved populations of massive stars. 

OO Key words. Stars: early-type - Stars: massive - binaries: spectroscopic - binaries: close - open clusters and associations: individual (30 Dor) - 
C*~) Magellanic Clouds 

<3\ 1. Introduction Type Ib/c supernovae (Yoon et al. 2010), short and long-duration 

gamma-ray bursts (Podsiadlowski et al. 2010). 

Massive binaries are spectacular systems that may have been 
^ among the first objects that formed in the early universe (Turk Among the various multiplicity properties, the most investi- 
> et al. 2009; Stacy et al. 2012). Because the binary fraction among gated one is the observed fraction of spectroscopic binaries, that 
massive stars is high (Mason et al. 2009) and as a large fraction sets a lower limit on the true binary fraction. Based on an exten- 
of them are part of short period systems (Sana & Evans 201 1), sive compilation of the literature, Mason et al. (2009) showed 
£d an accurate knowledge of their binary properties is crucial to un- that 5 1 % of the Galactic O-type stars investigated up to 2008 by 
derstand the role of massive stars as a population (Langer et al. multi-epoch spectroscopy are spectroscopic binaries. This frac- 
2008; de Mink et al. 2009). On the one hand, the fraction of tion is 56% for objects in clusters and OB associations. The spec- 
binaries and their orbital parameter distributions are the tracers troscopic survey of Galactic O and WN stars (Barba et al. 2010) 
of massive star formation and of the early dynamical evolution obtains a similar fraction: 60% of the 240 southern stars inves- 
of their host star cluster. On the other hand, they determine the tigated indeed display significant radial velocity (RV) variations 
frequency of evolved massive binary systems, including high- with an amplitude larger than 10 kms 1 . Finally, studies of in- 
mass X-ray and double compact binaries (Sadowski et al. 2008), dividual young open clusters, such as IC 1805 (De Becker et al. 

2006; Hillwig et al. 2006), IC 1848 (Hillwig et al. 2006), IC 

2944 (Sana et al. 2011), NGC 2244 (Mahy et al. 2009), NGC 

* Based on observations collected at the European Southern 6231 (Sana et al. 2008), NGC 661 1 (Sana et al. 2009) and Tr 

Observatory under program ID 182.D-0222 16 (Rauw et al. 2009) reveal observed binary fractions ranging 

** Hubble Fellow between 30 and 60%. These cluster to cluster variations are so 
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far compatible with random fluctuations expected from the size 
of the samples (Sana & Evans 201 1). 

The intrinsic period and mass-ratio distributions of the O- 
type binaries have remained poorly constrained until recently. In 
a pioneering work, Kobulnicky & Fryer (2007) attempted to con- 
strain these distributions using a sample of 900 RV observations 
of 32 O- and 88 B-type stars in Cyg OB2 and a Monte-Carlo 
method to correct for observational biases. Their solution was 
unfortunately degenerate leading the authors to fix the period 
distribution to an Opik law (i.e., a flat distribution in the loga- 
rithm of the separation or, equivalently, of the period). They find 
a very high intrinsic binary fraction of 80 to 100%, although the 
range of separations considered implies an extrapolation of their 
results over one to two orders of magnitude outside the sensi- 
tivity domain of their observations (see also Sect. 6.2). Based 
on a smaller sample of 13 O+OB and 37 B+B eclipsing bina- 
ries (Harries et al. 2003; Hilditch et al. 2005), Pinsonneault & 
Stanek (2006) suggested the presence of a population of equal- 
mass ('twin') binaries albeit Lucy (2006) argued that the consid- 
ered sample was too small to draw statistically significant con- 
clusions on the presence of a peak close to unity in the mass-ratio 
distribution. 

Sana et al. (2012) used a Monte-Carlo method to retrieve 
the intrinsic multiplicity properties of O-type stars in nearby 
Galactic open clusters. Based on as ample of 71 O-type ob- 
jects, they showed that the period distribution does not follow 
the widely used Opik law but is overabundant in short period 
systems and that the mass-ratio distribution is flat in the range 

0. 1 < M2IM1 < 1, ruling out both the presence of a twin popu- 
lation of equal mass binaries and random pairing from a classical 
initial mass function as a process to associate components in a 
binary. These results have strong implications for the evolution 
of massive stars as over two-thirds of the stars born as O star 
will, during their lifetime, interact with a nearby companion in a 
binary system. 

The clusters considered by Sana et al. (2012) have however 
relatively low masses, in the range 1,000-5,000 Mq. In view of 
the much more energetic environment provided by 10 s - 10 6 M Q 
clusters, containing up to thousands of massive stars, it is cur- 
rently unknown whether the binary fraction and the parameter 
distributions would remain unaffected. Addressing this question 
has implications in the context of super-star-clusters and cluster 
complexes observed well beyond the Local Group (e.g. Gieles 
etal. 2010). 

The Tarantula Nebula (or 30 Doradus) in the LMC is one of 
the largest concentrations of massive stars in the Local Group 
(Walborn & Blades 1997) with a total mass of -5.5 xlO 4 M Q 
(Andersen et al. 2009). Combined with its youth (the central 
cluster R136 is approximately 1 to 2 Myr; de Koter et al. 1998) 
and the presence of several sub-populations, indicative of vio- 
lent star forming activity in the past tens of millions of years, 
this complex system is our closest view of a starburst-like region 
in the Local Universe and an ideal laboratory to investigate the 
characteristics of massive stars. 

By targeting -800 massive stars the VLT-FLAMES 
Tarantula Survey (VFTS; Evans et al. 201 lb, Paper I) provides 
a near-complete census of the 30 Doradus region. The science 
drivers of this program are presented in e.g. Evans et al. (201 la) 
and de Koter et al. (2011). The study of the multiplicity prop- 
erties of this region is an important component of VFTS and to 
this end a multi-epoch observing strategy has been implemented 
to identify binaries with periods < 10 3 days (see Sect. 4.3), 

1. e. the systems in which the components are expected to inter- 
act through Roche-lobe overflow. Early serendipitous findings 
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Fig. 1. Top panel: Number of RV epochs obtained for each object. One 
Argus object, VFTS 1022, has been observed at 13 epochs and falls 
outside the plotted range. Middle panel: number of lines available for 
absolute RV measurements. Bottom panel: time elapsed between the 
first and last suitable epochs for RV measurement. The solid line corre- 
sponds to objects observed with Medusa while the dashed lines indicate 
objects observed with the Argus IFU. 



(Taylor et al. 201 1; Dufton et al. 201 1) clearly emphasize the im- 
portance of binaries and the large potential these systems have 
for our understanding of binary properties and evolution. 

In this paper, we use radial velocities (RVs, see Sect. 2) and 
a variability analysis to investigate the multiplicity properties of 
the 360 O-type stars in the VFTS sample. We first establish the 
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Fig. 2. Distribution of the O-type stars in the field of view. The popula- 
tion is dominated by the central cluster, NGC 2070, that contains 57% 
of our sample. NGC 2060, at 6.7' to the South- West, is slightly older 
and accounts for 22% of the sample, while the remaining 21% is spread 
throughout the field of view. The dashed circles, with a 2.4'-radius each, 
show the adopted regions for NGC 2070 and NGC 2060 (see Table 6). 



observed spectroscopic binary fraction in our sample (Sect. 3). 
Using a Monte-Carlo method, we correct for the observational 
biases and we constrain the intrinsic binary fraction and period 
and mass-ratio distributions (Sect. 4). Possible correlations of 
these quantities with spatial distribution, magnitude and spectral 
type are investigated in Sect. 5. Our findings are discussed in 
Sect. 6 and summarised in Sect. 7. 



2. Observations, sample and RV measurements 

2.1. Data overview 

The VFTS data and the data reduction process are extensively 
described in Paper I. Here, we provide a brief description of 
the VFTS campaign and observational settings related to the 
present study. Using nine different plate configurations, VFTS 
has collected multi-epoch fibre spectroscopy of about 800 early- 
type stars in the Tarantula region, covering a 25' diameter field 
of view centered on R136. Each observation consisted of two 
30-min exposures taken back to back. In this paper, we analyze 
the O-type star population observed with the Medusa fibres and 
the low resolution settings LR02 and LR03 of the Giraffe multi- 
object spectrograph. These data provide continuous coverage of 
the 3960-5070 A region at a spectral resolution of about 0.6 A. 
The higher resolution data obtained in the Ha region are not used 
in the present analysis. 

Most of the time sampling is provided by the LR02 obser- 
vations (3960-4565 A, A/AA « 7000), with five pointings cov- 
ering time scales of days to months and one additional point- 
ing taken about 300 days after the first observations. The LR03 
setup (4500-5070 A, A/AA = 8500) was observed three times 
with no specific time constraint. As a result some of the config- 
urations were observed consecutively to LR02 data while others 



are spread out and provide additional time sampling. For oper- 
ational reasons, one LR02 plate configuration - field C - was 
reobserved a seventh time in Oct 2010, extending the temporal 
baseline to almost 700 days for 48 O stars in our sample. The 
signal-to-noise ratio (S /N) is however lower than average and 
this last pointing provides little additional information. We have 
performed the complete analysis presented in Sect. 3 and 4 with 
and without these supplementary data and found no significant 
differences in the results. In the rest of the paper, we ignore that 
epoch in order to preserve the homogeneity of the temporal sam- 
pling for all stars in our data set. 

The combined effect of chromatic atmospheric diffraction 
and of seeing may strongly impact the fraction of the starlight 
that is injected into the Medusa fibres. As a consequence, the to- 
tal response of the observing system may present significantly 
different shapes from one exposure to another. Each 30-min 
spectrum is therefore normalised individually using the proce- 
dure described in Appendix A that provides a high degree of 
automatisation and guarantees an homogeneous handling of the 
data. 

Contiguous observations were then combined to improve the 
S /N, yielding a typical number of six epochs per targets and 
providing a minimum timebase of about 300 days for most of 
the stars. Statistics of the number of epochs, diagnostic lines for 
RV measurements and timebase are provided in Fig. 1. 

2.2. O-type stars in VFTS 

The sample analyzed in this work includes all O-type stars 
observed during the VFTS campaign and ignore the B-type 
stars, the Wolf-Rayet stars and the transitional objects to the 
Wolf-Rayet domain (the so-called "slash-stars"). The multiplic- 
ity properties of the B stars are investigated in a separate paper 
within the VFTS series (Dunstall et al., in preparation). 

The O-type stars were classified by comparison, at a resolv- 
ing power of 4000, with spectral standard stars. The detailed 
spectral classification is presented in Walborn et al. (in prepa- 
ration). In total, 332 O-type stars observed with Medusa form 
the bulk of our sample. 

The densest regions of 30 Dor, in and around the core cluster 
R136, were further observed with the Argus integral-field unit 
(IFU) and the LR02 Giraffe setup. The RVs of these objects 
were measured using a similar approach to the one described 
in Appendix B and the same set of lines and rest wavelengths 
(see Table B.l and Henault-Brunet et al. 2012). The Argus RV 
measurements of 28 additional O stars observed with the Argus 
IFU have been incorporated in the present analysis. In total, the 
sample that we analyze thus contains 360 O-type stars and is the 
largest homogeneous sample of O stars investigated for multi- 
plicity to date. 

Fig. 2 shows the spatial distribution of the O stars in the 
30 Dor field of view. About half of our sample (184 objects) is 
concentrated in a 2. 4' -radius around R136, the young and mas- 
sive cluster at the core of NGC 2070. Seventy-three targets are 
spatially associated to NGC 2060, a slightly older cluster at 6.7' 
South- West of R136 and the remaining 103 objects are spread 
across the field of view. 



2.3. Completeness 

The list of potential targets to be observed by the VFTS was se- 
lected based on a simple magnitude cut-off at V = 17 (see Paper 
I). The allocation of the MEDUSA fibres to VFTS sources was 
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Table 1. List of lines (marked with an 'x') used to obtain the RV measurements listed in Table 2. The full version of the table is available 
electronically. 
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Table 2. Journal of the observations, RV measurements and their lcr 
error-bars for the stars in our sample. The full version of the table is 
available electronically. 



VFTS 


HJD 
-2400000 


RV 

(kms- 1 ) 


0RV 

(kms- 1 ) 


014 


54815.791 


276.4 


1.6 


014 


54822.694 


289.9 


2.2 


014 


54859.720 


275.6 


2.4 


014 


54891.573 


274.7 


1.9 


014 


55113.818 


280.7 


2.4 


014 


54813.776 


279.8 


1.5 


016 


54817.734 


191.9 


0.7 


016 


54822.547 


187.6 


0.9 


016 


54860.594 


189.1 


0.9 


016 


54890.530 


188.7 


2.0 


016 


55112.850 


189.4 


1.2 


016 
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285.5 
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54827.772 


229.2 
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54860.643 


283.4 
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15.5 
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017 


54836.563 


260.2 


1.5 



subject to further constraints related to the physical size of the 
magnetic buttons supporting the fibres, precluding observation 
of all sources in the most crowded parts of the field-of-view. The 
VFTS completeness is about 70% of the total O-star population 
in the field of view and, importantly, does not show a bias to- 
wards the brighter end of the magnitude distribution. Beside the 
physical constraint on the fibre allocation, the VFTS sample is 



complete for the O-type stars, except for the late O-type objects 
that suffer from a significantly larger extinction - at least one 
magnitude - than typical in the field of view. 

2.4. RV measurements 

The radial velocities of our targets are obtained by measuring 
the Doppler displacement of spectral lines. Because the RVs are 
a key ingredient of this study, details of the method are provided 
in Appendix B, including alternative approaches that have been 
considered. Our method of choice, Gaussian fitting of selected 
sets of He lines, is prompted by four criteria: (i) robustness, (ii) 
internal consistency across the full range of O spectral sub-types, 
and the ability to provide (iii) absolute measurements and (iv) 
error bars. An extension to the standard line-by-line Gaussian 
fitting is obtained by simultaneously fitting all the data available 
for a given object (i.e., all the epochs and spectral lines). The 
different lines fitted at a given epoch are required to provide the 
same Doppler shift while the shape of the individual lines is kept 
constant through all the epochs. This approach provides two sig- 
nificant improvements compared to line-by-line fitting: a strong 
reduction of the number of parameters of the fit and a signifi- 
cant gain in robustness in the RVs measured from data with poor 
S/N. 

Special attention has been given to identify a set of lines that 
provide consistent RV measurements across the range of O spec- 
tral sub-types. Lines displaying systematic deviations as a result 
of, e.g., line blending or wind effects have been discarded (see 
Appendix B.2 for a detailed discussion). The best consistency 
is observed between RVs obtained from the Hen /L14200, 4541, 
He i /L14387, 4713 and 4922 lines. Whenever possible, we have 
thus obtained absolute RV measurements based on these five 
lines or a subset of these lines depending on the spectral type, 
wavelength coverage and degree of nebular contamination. The 
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Table 3. Overview of the results of the multiplicity analysis. Cols. 2 
and 3 provide the maximum significance of the RV variations and their 
maximum amplitudes, respectively. Col. 4 indicates whether LPV are 
detected. Col. 5 reports on the morphology of the lines and the last 
column indicates whether the targets are revealed as multiple sources 
by the HST The full version of the table is available electronically. 
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spectral lines used for each object are reported in Table 1 while 
Table 2 provides the journal of the observations and, for each 
epoch, the measured RV and lcr uncertainty. 

The typical number of spectral lines used in measuring the 
RVs is between two and five (Fig. 1). For VFTS 051, 400, 
444, 453, 531, 565 and 587, none of these lines were of suf- 
ficient quality. RVs for these targets were measured using the 
He i+ii /14026 and/or He n /14686 lines. Because these are late O- 
type stars the measured RVs are probably only weakly affected 
by temperature and/or wind effects but caution should be applied 
in using the measured RVs in absolute terms. A number of ob- 
jects displayed double-lined profiles (see Sect. 3.2.1). Measuring 
RVs for both components sometimes required the use of lines 
that are not suited to provide absolute RVs. Because SB2s are 
not suited for absolute RVs unless a full orbit determination can 
be achieved, we relax our constraints for these objects and in- 
clude other He lines as indicated in Table 1 . For the SB2 cases, 
Table 2 lists the RV of the earliest component of the pair. 

3. Multiplicity analysis 

In this section we investigate the multiplicity of the stars in our 
sample using various approaches that combine variability anal- 
ysis (Sect. 3.1) and the search for composite sources (Sect. 3.2) 
revealed either by multiple signatures in their spectrum or by 
HST imaging. Sect. 3.3 establishes the observed spectroscopic 
binary fraction and Sect. 3.4 compares our results with earlier 
studies. The results of the multiplicity analysis are summarised 
in Table 3, the full version of which is available online. 



Table 4. Overview of the results of various variability and multiplicity 
criteria used in this paper 



Variability 


Constant 


Weak variability 


Variable 


RV variability 


195 


39 


126 


TVS analysis 


253 


15 


96 


Multiple signature 


Single/Isolated 


Candidate 


Multiple 


Profile morphology 


261 


24 


63 


HST images" 


303 


5 


24 



NOTES: a. Nineteen stars have no HST images. 

3.1. Variable sources 
3.1.1. RV variability 

To identify objects displaying significant RV variations, we use a 
statistical test. The null hypothesis of constant RV is rejected if, 
for a given star, any two RV measurements deviate significantly 
from one another, i.e.: 



0"detect = max 



\Vi - Vj 



2 + CT 2 
J J 



>4.0, 



(1) 



where v, and <x, are the RVs and their lcr errors for a given object 
at epoch ;'. The confidence threshold of 4.0 provides a probabil- 
ity of about one false variability detection in every 1000 objects, 
given the sampling and accuracy of our measurements. This cor- 
responds to an expected impact on the measured fraction of RV 
variable objects of 10~ 3 . We find that 165 sources (46%) display 
significant RV variability. 

Keeping all other things equal, Eq. 1 provides a more sensi- 
tive criterion to detect variability than a comparison of the v, 's 
with the average RV of the object. The results of Eq. 1 are also 
equivalent to those of a variability test based on the goodness of 
fit of a constant RV model (see e.g. Henault-Brunet et al. 2012) 
for over 96% of our sample. 

3.1.2. Line profile variability 

We also search for line profile variability (LPV) using a Time 
Variance Spectrum (TVS) analysis (Fullerton et al. 1996). The 
TVS analysis is sensitive to variations in the amplitude, shape 
and position of the line profile, as well as to inconsistency in 
the nebular correction between the various epochs and to techni- 
cal defects such as cosmic rays and continuum localisation. The 
method also allows us to verify the consistency of the normalisa- 
tion procedure described in Appendix A. We adapt the original 
TVS approach of Fullerton et al. (1996) to take into account the 
known error bars at each pixel, rather than to rely on the overall 
S IN of the spectrum. This is needed because the S /N is variable 
as a function of wavelength and depends on the chosen setup. 
We thus define : 



TVS (A) = 



N 



N ■ 



^ (MX)- < F(A) >y 



1 



(2) 



where Fj and <x, are the observed normalised flux and its lcr error 
as a function of wavelength A at epoch i. < F > is the weighted 
average flux computed over all available epochs. The statistical 
significance threshold is adopted to be a = 0.01, i.e.: 



TVS (A) > 



N 



N — 1 



1}/ ?^)' 



(3) 
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where P x i{a,N - 1) gives the cutoff value in a x 1 distribution 
with N - I degrees of freedom such that the probability that 
a random variable exceeds P x i{a,N - 1) is equal to a. In our 
data, most of the variability detected through the TVS analysis is 
produced by non-optimum nebular correction. It results from the 
combined effects of the limited spectral resolution, the spatially 
variable level of nebular emission and the fact that FLAMES 
does not provide a nebular spectrum close to the target. 

In addition to residuals of the nebular correction, the TVS 
detects significant variability in 96 cases. By far, most cases (92) 
are also associated with significant RV variations so that the TVS 
analysis does not significantly add to the number of detected 
variable stars. 

3.2. Composite sources 

3.2.1. Double-lined binaries 

The spectra obtained for each target were visually inspected to 
search for the signature of double-lined and asymmetric profiles. 
Despite the limited S /N of the individual epoch spectra, double- 
lined profiles were clearly identified in 50 cases. Higher multi- 
plicity profiles were not observed. Asymmetric or variable pro- 
files are observed for 13 objects. Candidate SB2s and objects for 
which asymmetric profiles are suspected represent another 21 
cases. For pronounced SB2s, RV measurements were attempted 
using two Gaussians per profiles. Reliable fits could be achieved 
in 43 of the 50 cases. In two cases the stronger component - 
likely the primary, more massive star - shows no RV variation. 
These objects are thus added to the list of binary candidates al- 
though one cannot confirm the physical link between the two 
component of the pairs. 

3.2.2. HST observations 

We investigate the targets for visual companions using HST data 
taken with WFC3/UVIS as the primary camera with ACS/WFC 
in parallel, both using the F775W filter. These data constitute 
the first epoch of a proper motion program in the region (GO 
12499, PI: Lennon, see de Mink et al. 2012; Sabbi et al. 2012). 
It provides a ~ 16' x 12' mosaic centered on R136 covering 
over 90% of the O stars in the Tarantula Survey. We find that 19 
O stars observed with Medusa and six with ARGUS have one 
or several companions that are close and bright enough, for the 
FLAMES spectra to be a composite of multiple sources. These 
and an additional five multiple candidates are listed in Table 3. 

The HST point spread function, of ~0.1\ corresponds to 
a physical separation of about 5 x 10 3 AU at the distance of 
30 Dor. HST thus probes a very different, non-overlapping, sep- 
aration range compared to the spectroscopic binaries discussed 
in Sect. 3.1. Five of the visual pairs identified by HST also show 
significant RV variations, indicating that at least one of the com- 
ponents of the visual system is also a spectroscopic binary 

3.3. Observed binary fraction 

Of all the diagnostics discussed above, RV variability is by far 
the most robust test to identify variable sources in our data. 
Complementary tests (Sects. 3.1.2 to 3.2.2) do not significantly 
change the number of detected systems. Furthermore the RV 
variability tests defined by Eq. 1 do not rely on subjective ap- 
preciations such as e.g., the visual inspection of the line profile 
morphology. This is a necessary condition to model the observa- 
tional biases by means of Monte-Carlo simulations (Sect. 4.1). 



0.5 F 




0.1 r .■""[" i 

0.0 E - .................. : 

20 40 60 80 100 

C fkm s~ ) 

Fig. 3. Variations of the fraction of systems below and above the crit- 
ical RV variation amplitude C separating the spectroscopic binaries 
(solid line) from low-amplitude RV variables (dotted line). The vertical 
dashed-dotted line indicates the adopted threshold of C = 20 km s -1 . 



In the following we will thus focus on spectroscopic systems that 
reveal themselves through their RV variability. 

The test provided by Eq. 1 only identifies significant RV 
variations. To distinguish cases where the observed variability 
is genuinely caused by orbital motion in a binary system from 
potential intrinsic variability due to e.g. photospheric or wind 
effects, we impose a minimum amplitude threshold C on any sig- 
nificant deviations computed from any individual pair of epochs. 
I.e., an object is considered a spectroscopic binary if at least one 
pair of RV measurements satisfies simultaneously 

|v; - V/| 

J > 4.0 and \v t - vj\ > C. (4) 



Fig. 3 shows the influence of C on the identified binary frac- 
tion. It reveals a kink at about 15 kms _I that may mark the 
transition between a regime were photospheric variability has 
a significant contribution (ARV < 15 kms~') and the regime 
of genuine spectroscopic binaries. A crude extrapolation of the 
contribution of the binaries into the low-amplitude ARV regime 
suggests the presence of about 5% of genuine binaries with 
ARV < 15-20 kms~'. Similarly, we estimate that photospheric 
variabilities affect as little as 6 to 7% of our sample. 

In Section. 4 we will correct for the observational biases and 
consider the fraction of long period systems that potentially fall 
below the ARV cutoff. As a consequence, the exact value of the 
adopted cutoff is of little importance but has to be large enough 
to avoid a significant number of false detections due to measure- 
ment errors and/or photospheric variability. 

Because photospheric variations in supergiants can mimic 
variations with amplitudes of up to 20 kms~' (e.g. Ritchie et al. 
2009), we conservatively adopt C = 20 kms~'. Eleven per cent 
of the sample (39 stars) are thus identified as variable but fea- 
ture an amplitude of the significant RV variations that remains 
below 20 km s" 1 . These objects might be long period binaries or 
stars with photospheric variability and, as discussed above, we 
estimate an equal contribution of both categories. 

Given Eq. 4 with C — 20 kms~', the observed spectroscopic 
binary fraction identified through RV variation of the component 
with the strongest line - likely the primary - is thus 35.0 ± 2.5%, 
where the lcr error bar gives the statistical error that results from 
the size of the sample (Sana et al. 2009). 
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Fig. 4. Histogram of the peak-to-peak amplitude RV variations of each 
object. Only significant variations according to Eq. 4 are considered. 
Therefore the first bin is empty. 
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Fig. 5. Cumulative distribution of the minimum time difference AHJD 
for a binary to display significant RV variations according to Eq. 4 and 
for different RV amplitude thresholds (C) in the RV variation ampli- 
tudes. The cumulative distributions are normalised to the total number 
of binaries. 



3.4. Comparison with earlier results 

Bosch et al. (2009, BTT09) reported on a similar study of 54 O 
and early B stars in NGC 2070. They obtained R ~ 3700 seven- 
epoch spectroscopy covering a wavelength range equivalent to 
the coverage of our LR02+LR03 settings. Using the external- 
to-internal velocity dispersion ratio (TeIci and a threshold of 3 
they identified 17 RV variables and nine additional SB2 bina- 
ries and binary candidates for a maximum detected binary frac- 
tion of 0.48. The results of BTT09 thus suggest a larger binary 
fraction than the one we observed in VFTS. To understand this 
difference, we perform a detailed comparison of the VFTS and 
BTT09 approaches. 

We first apply our variability and binary criteria on BTT09's 
RV measurements. We detect 24 RV variable objects among the 
48 objects listed in their Table 2. The difference comes from 
the fact that the ratio cr^lcr^ is artificially reduced if one mea- 
sure has a significantly worse accuracy that the others. All of the 
stars identified as RV variables have peak-to-peak amplitudes in 
excess of 20 kms -1 and thus qualify as spectroscopic binaries 
according to the criteria of Eq. 4. Applying our method on the 
BTT09 RV data set, we measure a spectroscopic binary fraction 



Table 5. Overview of properties of the Monte-Carlo grid. From left to 
right we list the physical parameter, its corresponding probability den- 
sity function (pdf) and domain of applicability, the name of the pdf vari- 
able and its allowed range and step size. The last row gives the allowed 
range and step size of the binary fraction. 



Parameter 


pdf 


Domain 


Var. 


Range 


Step 


P 


(logjoW 


0.15-3.5 


n 


-2.50- +2.50 


0.05 


1 


<t 


0.1 - 1.0 


K 


-2.50- +2.50 


0.05 


e 


e" 


10- 5 - 0.9 


1 


-0.5 (fixed) 


n/a 


/bin 


n/a 


n/a 


/bin 


+0.20- +1.00 


0.01 



of 50.0% while BTT09 obtained 35.4% by using the ^/cn ra- 
tio. Two of the seven binaries not identified as RV variable by 
BTT09 are still accounted for in their total binary count because 
of asymmetry in their line profiles. Even though we require a 
stricter (Act) threshold, we conclude that our binary criteria are 
slightly more sensitive than the one used by BTT09 based on 
cr E /cri > 3. 

As a second step in comparing our results with those of 
BTT09 we compute the binary fraction from the overlapping 
sample between BTT09 and the present work, which are 34 of 
the 54 objects studied by BTT09. Of these 34 objects, BTT09 
reported RVs for 30 of them among which they identified 12 
binaries (40%) through RV variations. Among the same set 
of 30 stars, VFTS identify 13 binaries (43%) and 6 (20%) 
low-amplitude variable objects, with ARV ranging from 7 to 
19kms-'. 

On the overlapping sample, the two studies agree well, with 
VFTS being slightly more sensitive due to the improved RV 
accuracy and statistical criteria. The significantly lower binary 
fraction observed in the whole VFTS sample relative to the work 
of BTT09 is attributed to the focus of BTT09 on the brightest 
stars. We will show in Sect. 5.2 that fainter O-type objects have 
an intrinsically lower binary fraction. 



4. Intrinsic multiplicity properties 

The detected binary population provides a minimum estimate of 
the true binary population in the Tarantula region. The number of 
undetected binaries depends on the sensitivity of our campaign 
in terms of RV accuracy and time sampling and on the intrinsic 
distribution of the orbital parameters. The problem of correcting 
for the observational biases is thus ill-posed and we employ here 
a Monte-Carlo approach to estimate the impact of the observa- 
tional biases and to provide constraints on the intrinsic binary 
fraction and distributions of orbital parameters. 

4. 1 . Monte-Carlo method 

In our Monte-Carlo approach, we simulate massive star popula- 
tions with intrinsic multiplicity properties randomly drawn from 
adopted parent distributions. Accounting for observational bi- 
ases due to sampling and measurement uncertainties we search 
for sets of distributions that reproduce the properties of the ob- 
servations in three aspects: the observed binary fraction, the 
peak-to-peak amplitude ARV of the RV variations (Fig. 4) and 
the minimum time scale AHJD for significant RV variation to be 
observed (Fig. 5). 

The method is detailed in Appendix C. In short, simu- 
lated and observed distributions are compared by means of 
Kolmogorov-Smirnov (KS) tests. The detected binary fraction 
predicted by the simulations is compared to the observed frac- 
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Fig. 6. Projections of the global merit function H' on the various pairs of planes defined by n, k and / bin . The absolute maximum is indicated by a 
cross (+). Contours indicate loci of equal-values of the merit function, expressed as a fraction of its absolute maximum (see legend). 



tion using a Binomial distribution (Sana et al. 2012). A global 
merit function (H') is constructed by multiplying the two KS 
probabilities obtained for the ARV and AHJD distributions and 
the Binomial probability: 

H' = Pks(MV) x Pks(AHJD) x B(N hin , N, f£™ 1 ), (5) 

where / b s ™ ul is the simulated fraction of detected binaries, Nu n 
is the detected number of binaries in the observations and N is 
the population size. 

Following Kobulnicky & Fryer (2007) and Sana et al. (2012), 
we use power-laws to describe the intrinsic distribution func- 
tions of orbital parameters: /(log 10 P/d) ~ (logjoP)^, f(q) ~ q" 
and /(e) ~ e n . The exponents n, k and 77 are then varied to ex- 
plore different shapes for the parameter distributions. Because 
the global merit function is insensitive to the eccentricity distri- 
bution (see Fig. C. 1), we cannot constrain 77 and we adopt the ec- 
centricity distribution measured using the O star populations of 
Galactic young open clusters (Sana et al. 2012), i.e. /(e) ~ e~ 05 . 
We thus explore the behaviour of H' in the three-dimensional 
parameter space defined by the intrinsic binary fraction and dis- 
tributions of the orbital periods and mass-ratios. The exploration 
domains and mesh sizes are given in Table 5. 

4.2. Results 

Fig. 6 displays the merit function (H') projected on planes de- 
fined by the various pairs of parameters. It shows that the opti- 
mum is well defined within the explored ranges. The best agree- 
ment between the observed and simulated distributions is ob- 
tained for n = -0.45 + 0.30, k = -1.00 ± 0.40 and a intrinsic bi- 
nary fraction of /,;„ = 0.5 1 ± 0.04. The quoted errors have been 
estimated as the lcr confidence interval on the retrieved parame- 
ters of 50 synthetic data sets randomly drawn from parent distri- 
butions close to our best fit distributions (see Appendix C.3 and 
the lower part of Table C.l). The comparison between the sim- 
ulated and observed cumulative distribution functions (CDFs) 
of the peak-to-peak RV amplitudes and of the variability time 
scales is presented in Fig. 7. 

While the overall shapes of the observed and simulated func- 
tions are well matched, some features of the observed distribu- 
tions are not fully reproduced by the simulations. The CDF of 
the peak-to-peak RV amplitudes presents an overabundance of 
systems with ARV in the range 300 to 400 km s . This could be 
reproduced by increasing the fraction of short period systems, 
i.e. a more negative n exponent, but the CDF of the variability 
time scale would then not be properly reproduced. Alternatively, 
a flatter distribution of mass-ratios would help to improve the 



CDF(A/?V) fit in the 300 to 400 km s -1 range but would overes- 
timate it in the lower amplitude range (ARV < 100 kms~'). 

We also attempt to reproduce the overabundance of systems 
with ARV in the range 300 to 400 kms~' by including specific 
populations of short period systems, of equal mass systems, of 
high-mass primaries but none of these, nor any combinations of 
them, were able to reproduce the observed localised bump. 

The bump in the CDF(A7?V) curve comprises less than 15 
systems and can also result from statistical fluctuations. Indeed, 
the individual Pks probabilities are already large: Pks(ARV) 
0.4 and P K s(AHJD) * 0.8. We thus consider the present fit to 
provide an adequate representation of the data. 

The value of the period distribution index is remarkably 
close to the value 7Tgoc = -0.55 ± 0.22 found in the O star pop- 
ulation of Galactic open clusters (GOCs, Sana et al. 2012). The 
intrinsic binary fraction is 51% in VFTS compared to 69 ± 9% 
in the GOCs. Finally, the mass-ratio distribution obtained in 
the present study favours lower mass companions while the 
mass ratios in the Sana et al. study were equally distributed 
(kgoc = -0.1 + 0.6). Both studies do however agree within 2cr. 
The implications of these results are discussed in Sect. 6. 



4.3. Sensitivity domain of the VFTS campaign 

We use the method of Sana et al. (2009) to estimate, given the 
distributions of the orbital parameters that we find, the sensitivity 
of our campaign with respect to the orbital period, the mass ratio, 
the eccentricity and the primary mass. 

Compared to the Sana et al. (2009) approach, we now in- 
clude the measurement uncertainties obtained at each epoch in 
the simulations and we implement the binary detection criteria 
given by Eqs. 1 and 4 to consistently follow the binary detection 
strategy laid out in this paper. 

The detection probability curves in Fig. 8 are projections 
from the multi-dimensional space onto each of the four dimen- 
sions investigated. The binary detection probability is a strong 
function of orbital period. It moderately depends on the mass 
ratio except for systems with a large mass difference. The de- 
tection probability is also weakly affected by the primary mass 
as well as by the eccentricity. Because we have not taken into 
consideration the fact that high eccentricity systems may prefer- 
entially have a long period, we underestimate the detection rate 
at very high eccentricities (e > 0.7). 
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Fig. 7. Comparison between the observed (crosses) and simulated (lines) cumulative distributions of the peak-to-peak RV amplitudes (top row) 
and of the variability time scales (bottom row). In the left- (resp. right-) hand panels, we vary the exponent n of the period distribution (resp. k of 
the mass-ratio distribution) by ±lcr. The upper-left legends indicate the values of n, k and 77 considered in each panel. The bottom-right legends 
give the overall VFTS detection probability for the adopted parent distributions. 



5. Binary fraction 

In this section, we discuss how the binary fraction correlates with 
the spatial distribution, brightness, spectral-type and luminosity 
class of the stars in our sample. Supplementary figures support- 
ing this discussion are provided in Appendix D. 



5.1. Spatial variations 

We search for variations in the observed binary fraction through- 
out the field of view using various techniques. First we divide the 
stars in three samples based on their spatial location: stars within 
a 2.4' radius from the centre of NGC 2060 and NGC 2070, 
and stars outside these two clusters. NGC 2060 tends to dis- 
play a slightly larger fraction of low-amplitude RV variable stars 
(Table 6). This suggests a larger number of objects affected by 
photospheric variability, likely supergiants, reflecting the older 
age of the cluster, hence the larger fraction of supergiants (see 
also Fig. 13). When the sample is limited to stars with large am- 



Table 6. Observed fraction of RV variables and binaries among different 
O star populations in 30 Dor. 



Population (size) RV variables (Eq. 1) Binaries (Eq. 4) 
C(kms-') 0.0 20.0 



Overall (360) 


0.458 


± 0.026 


0.350 


± 0.025 


NGC 2070° (184) 


0.429 


± 0.036 


0.332 


± 0.035 


NGC 2060* (73) 


0.562 


± 0.058 


0.384 


± 0.057 


Remaining 1 ' (103) 


0.437 


± 0.049 


0.359 


± 0.047 



notes: a. Denned as a 2.4' radius around R136 (a = 05 : 38 : 42.0, 
6 = -69 : 06 : 02); b. Defined as a 2.4' radius around 
a = 05 : 37 : 45.0, 6 = -69 : 10 : 15; c. All stars not included in 
NGC 2060 nor NGC 2070. 



plitude RV variations, the genuine binaries, there are no obvious 
spatial differences between the various populations. 

We also look for spatial variations in the field of view us- 
ing a 2D smoothing. We do not identify significant variations 
save for two trends. There are slightly fewer binaries towards 
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Fig. 10. From left to right, distributions of the RV constant, likely single stars, the binary candidates (low-amplitude RV variables) and the binary 
stars in the 30 Dor field of view. 
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Fig. 8. Overview of the binary detection probability achieved by the 
VFTS campaign as a function of the orbital period, the mass ratio, the 
eccentricity and the primary mass. [See Sect. 4.3 for more details.] 



the Southern and, to a lesser extent, Northern parts of the field 
of view, as well as some enhancement in between NGC 2070 
and NGC 2060. These trends are significant at the 95% level 
(« l.lcf), but do not reach the 2cr confidence level. 

Finally, we also build the radial cumulative distribution of 
the binary fraction, starting both from the inside and from the 
outside (Fig. 9). The centre of the field of view is taken at the 
core of R136, as defined in Table 6. Overall, there is no signifi- 
cant variability but for the edges of the field. The binary fraction 
for the 22 outermost objects is 0.14 ± 0.07 which is just at 3<x 
from the average binary fraction in the whole sample. A visual 
impression of the lack of binaries in the outer regions can be ob- 
tained by comparing the locations of the single and binary stars 
as shown in Fig. 10. 

5.2. Brightness dependence 

Objects in our sample span five magnitudes in the V band. We 
use the V-band photometry published in Paper I to investigate 
any variations of the binary fraction with magnitude. We also 
employ /T s -band photometry from the VISTA Magellanic Clouds 
(VMC) Survey (Cioni et al. 201 1) to the same end. Specifically, 
we use PSF-fitting photometry of the 30 Dor VMC Survey tile 
from Rubele et al. (2012), with cross-matches to the VFTS 
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Fig. 9. Cumulative radial distribution of the observed binary fraction 
as a function of the distance to the centre of the field. The black solid 
line shows the distribution computed inside-out, and the grey/red line 
shows the distribution computed outside-in. The dashed curves indicate 
the ±1<t confidence interval on both curves. 

sources listed by Zaggia et al. (in preparation). A^-band magni- 
tudes from the PSF-fitting analysis are unavailable for 15 VFTS 
stars, either due to crowding, blending or saturation issues. VMC 
photometry is also missing for the Argus targets and these are 
henceforth excluded from the present analysis. Where possible, 
the VMC Survey results are supplemented with 2MASS pho- 
tometry using the cross-matches given in Paper I. 

The analysis of both bands provides similar results. Because 
of the variable extinction in the region (Mafz Apellaniz et al., 
in preparation) we focus here on the K s data only. Fig. 1 1 dis- 
plays the magnitude distributions of the VFTS O stars and the 
observed binary fraction in various magnitude bins. The lat- 
ter seems relatively uniform for objects brighter than K s = 
15.5 mag. For dimmer objects, the observed binary fraction 
starts to decrease and reaches 16% at the magnitude cut-off of 
our survey: K s ss 16.5 mag, corresponding to V = 17 mag. 

The apparent lower binary fraction for fainter stars can be 
explained by the combination of two observational effects: 

(i) spectra from fainter stars have a lower S /N, thus RV mea- 
surements are on average less accurate and small RV varia- 
tions are more difficult to detect; 
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Fig. 11. Top panel: Number of objects (dashed line) and of binaries 
(solid line) per K s -band magnitude bin. Bottom panel: Observed binary 
fraction in corresponding magnitude bins (solid line) and intrinsic bi- 
nary fraction (dashed line). Error bars are plotted whenever a bin con- 
tains more than one star. Fig. D.2 provides similar figures showing the 
dependence with luminosity class. 

Table 7. Number of objects, observed binary fraction detection proba- 
bility and intrinsic binary fraction as a function of the A" s magnitude. 



K s mag range 


Size 


xobs 
J bin 


^detect 


/bin 


< 14.0 


50 


0.345 ±0.106 


0.731 


0.472 


14.0-14.5 


34 


0.405 ±0.129 


0.705 


0.575 


14.5-15.0 


42 


0.449 ± 0.106 


0.676 


0.664 


15.0-15.5 


68 


0.449 ± 0.084 


0.635 


0.706 


15.5-16.0 


53 


0.245 ±0.115 


0.609 


0.403 


16.0-17.0 


24 


0.148 ±0.177 


0.572 


0.259 



(ii) binaries are typically brighter than single stars of the same 
spectral type by up to 0.75 mag, depending on the bright- 
ness of the companion. Binary systems thus tend to avoid 
the faintest bins, which then appear to present a lower binary 
fraction. 



Table 8. Observed fractions of low-amplitude RV variables and binaries 
as a function of the luminosity class. 



Luminosity class (size) 


low-ampl. RV variables 


Binaries 


V to I (293) 


0.11 


0.38 


V/TV (186) 


0.06 


0.38 


III (77) 


0.13 


0.42 


II/I (30) 


0.40 


0.23 
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Fig. 12. Top panel: Number of objects (dashed line) and of binaries 
(solid line) per O spectral sub-type. Bottom panel: Observed binary 
fraction as a function of the spectral sub-types. Error bars are plotted 
whenever a bin contains more than one star. Fig. D.3 provides similar 
figures for the different luminosity classes. 



To investigate the effect of (i) and estimate the intrinsic bi- 
nary fraction in each bin, we compute the average detection 
rate achieved in different magnitude bins using the Monte-Carlo 
method of Sana et al. (2009). We adopt the orbital parameter 
distributions obtained previously (Sect. 4.1). However, we do 
not take into account the likely dependence of the primary mass 
with the magnitude, i.e. we keep the primary mass range fixed 
to 15-80 Mq for all considered intervals. This impacts by a few 
per cent the detection probabilities listed in Table 7. When the 
results are computed for all of the O stars, a slightly different 
binary fraction is found than obtained in Sect. 4.1. Still, Table 7 
indicates that the lack of binaries at K s > 15.5 is an intrinsic 
effect, which can not be explained by the lower probability to 
detect binaries among the fainter objects in our sample. 

To investigate the effect of (ii) we simulate the magnitude 
distribution of an O-type population made-up of 50% single stars 
and 50% binaries. We find that the intrinsic binary fraction in the 
brightest bins was typically close to 60% but drops to 20-25% in 
the last bin, in excellent agreement with the intrinsic binary frac- 
tion measured in Table 7. Adding a slightly variable extinction 



easily allows us to extend the effect over the last two bins. If 
early B-type stars with a similar binary fraction are included in 
the simulations the magnitude bins at K s -45.5-16. 5 do not show 
such a drop. 

We thus conclude that the variation of the intrinsic binary 
fraction with magnitude is simply a consequence of the fact that 
our sample is defined by a cut-off in spectral type, i.e. the O- 
type sub-sample of the VFTS is not magnitude limited, save for 
the few stars with severe extinction. We also conclude that O 
star binary fractions obtained from magnitude-limited samples 
can be significantly overestimated, as such a selection criterion 
picks up extra binaries near the magnitude limit. 

This depletion of binaries for fainter objects is actually re- 
sponsible for the apparent difference between our overall results 
and the observed binary fraction obtained by Bosch et al.. Our 
observed binary fraction among objects with K s < 15.5 is 41 .4%. 
Based on the results of Table 7, we estimate an intrinsic binary 
fraction of 61.4% in that brightness range. 
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Fig. 13. Spatial distribution and multiplicity properties of the super- 
giants in the field of view. 



5.3. Spectral sub-type and luminosity class dependence 

To investigate possible variations of the detected binary fraction 
with the spectral sub-type and luminosity class, we ignore all ob- 
jects which spectra have a too poor S /N or are too heavily con- 
taminated by nebular lines to yield decent spectral classification. 
We also ignore the Argus sources that have a more limited spec- 
tral coverage and thus a less accurate spectral classification, and 
the seven Onfp stars, that present peculiar spectra as noted by 
the qualifier 'p' appended to their spectral type. We are left with 
295 out of 328 Medusa objects (i.e., 90% of the Medusa sam- 
ple). These objects are distributed among the luminosity classes 
as follows: 30 I-II objects, 77 III objects and 186 IV-V objects. 



Fig. 12 shows the observed binary fraction as a function of spec- 
tral sub-type and Fig. D.3 splits Fig. 12 according to luminosity 
class. Table 8 indicates the overall low-amplitude RV variable 
and binary fraction in each category. Most of the fluctuations 
are compatible with the sample sizes so that the binary fraction 
seems homogeneous across the spectral sub-types in the various 
samples. With two binaries and one low-amplitude RV variable 
object in a sample of seven, the binary fraction of the Onfp stars 
is also compatible with that of the main sample. 

The 06 stars display a significantly lower binary fraction 
than the average of the population. Inspection of Fig. D.3 reveals 
that the drop in the binary fraction at 06 is due to the dwarfs, 
with only one out of 1 3 06 V stars being detected as a binary, 
thus a binary fraction of 0.08 ± 0.07. 

The binary fraction among supergiants (23%) is smaller than 
among dwarfs and giants (~40%, Table 8). Interestingly, these 
numbers are roughly in line with expectations from binary evo- 
lution. The larger radii of I and II stars may imply that the short- 
est period binaries (P < 2-3 d) have already interacted. Because 
post-interaction systems are most likely to appear as single stars 
in RV studies (de Mink et al., in preparation), the fraction of 
detected binaries is thus expected to be lower. 

With only one binary out of 12 objects (corresponding to an 
observed binary fraction of 0.08±0.08) but a large fraction (42%) 
of objects showing low-amplitude RV variability, the 09.7 I and 
II stars are actually responsible for most of the differences be- 
tween the binary fraction among supergiants and among dwarfs 
and giants. Fig. 13 compares the spatial location of the 09.7 I-II 
objects with the other class I-II objects in the field. It reveals that 
the 09.7 supergiants preferentially belong to NGC 2060 while 
the hotter supergiants are distributed according to the ratio of the 
populations of NGC 2060 and NGC 2070. The data associated 
with the 09.7 supergiants are of a high quality, therefore exclude 
the possibility that the sample suffers from a higher detection 
threshold than the other supergiants. It may thus be interesting 
to speculate that the 09.7 I-II population in NGC 2060 contains 
a large fraction of apparently single post-interaction stars. 

5.4. Summary 

The binary fraction displays a high degree of homogeneity 
across the different populations. Only three sub-categories de- 
viate significantly from the mean: 

i. The overall binary fraction is lower in the outer region (r > 
7.8') of the field of view ; 

ii. The binary fraction is lower for the fainter stars in the sam- 
ple (K s > 15.5), an apparent effect that results from obser- 
vational biases and the definition of the sample in terms of 
spectral type; 

iii. The 09.7 I-II stars, mostly located in NGC 2060, are pre- 
dominantly single, though a relatively large fraction show 
low-amplitude RV variability. 

6. Discussion 

6.1. Nature of the pairing process 

We investigate two different pairing processes to associate the 
components of the binary systems and we compare the results 
of the mass-ratio distributions of simulated populations with the 
one derived from the observations. 

i. Random pairing: Primary and secondary masses are inde- 
pendently drawn, between 0.05 and 80 M G , from a mass 
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q = M 2 /M, 

Fig. 14. Mass-ratio distributions obtained using different pairing mech- 
anisms (solid and dashed lines) compared to those obtained with a 
power-law density distribution with different k values (dotted lines). 

function with a slope of -2.35. Stars are randomly associ- 
ated in pairs; 

ii. Truncated secondary mass function: The primaries masses 
are drawn as above. The secondary masses are drawn from a 
mass function truncated at the primary mass. 

To mimic an O star selection criterion, only pairs for which at 
least one of the components has a mass larger than 15 M Q are 
kept. In case (i), the positions of the primary and secondary are 
exchanged if the latter is the most massive component of the pair, 
so that q < 1 is respected. The q distribution is then truncated at 
q — 0.1 to confine to the same parameter range as discussed in 
this paper. 

Fig. 14 illustrates the mass-ratio distributions obtained. None 
of the tested pairing processes produces a distribution function 
that resembles f(q) oc q K with -1 < k < 0. It also shows that, 
as expected, random pairing results in k — -2.35 and that the 
truncated secondary mass function comes close to k — 1.0. 

6.2. On the universality of the multiplicity properties 

Our results can be discussed in the broader context laid out by 
recent studies of two different O star samples in the Milky Way. 
First, Kobulnicky & Fryer (2007) attempted to model the mul- 
tiplicity properties in Cygnus OB2 using the results of a multi 
epoch RV campaign targeting 120 OB stars (900 epochs in to- 
tal). Because they only adjusted the x~ distribution, the authors 
obtained a degenerate solution in the exponent of the period and 
mass-ratio distributions and they chose to fix ji at 0.0. The sep- 
aration regime considered in their modelling extends at least up 
to 100 A.U. and, for some of their simulations, to 10000 AU, 
i.e, quite beyond the sensitivity range of their spectroscopic ob- 
servations. In the considered parameter ranges, Kobulnicky & 
Fryer obtained a large binary fraction, at least 80%, and k values 
ranging from -0.6 to 1.0. Because of the different approaches 
and the different parameter ranges considered, direct compari- 
son with the current results is hazardous. 

Second, Sana et al. (2012) studied the multiplicity properties 
of the O star populations in nearby young open clusters. They 
performed a similar Monte-Carlo analysis as done here. Most of 
the binaries identified in their sample have constraints on their 
orbital parameters. This allowed them to test directly the period 
and mass-ratio distributions. They obtained an intrinsic binary 



fraction of 0.69 ± 0.09 and values of n = -0.55 ± 0.22 and k = 
-0. 1 ± 0.6 for the power-law exponents of the period and mass- 
ratio distributions. 

The period distributions obtained in both studies are in ex- 
cellent agreement. The mass-ratio distribution in Sana et al. is 
almost fiat while that of the VFTS sample is more abundant 
in lower mass-ratios. Both results are not inconsistent within 
2cr and, as discussed in Sect. C.2, the present approach is only 
weakly sensitive to k. The main difference resides in the intrinsic 
binary fraction which seems significantly lower in the Tarantula 
region. This aspect will be discussed in Sect. 6.3. 

6.3. Evolutionary implications 

Assuming that today's distributions are representative of the dis- 
tributions at birth and using a series of assumptions on critical 
period and mass-ratio values that delimit various binary inter- 
action scenarios, the integration of the distribution functions of 
Sect. 4.1 provides an estimate the evolutionary fate of the 30 Dor 
O star population. For comparison purposes, we adopt the same 
set of assumptions as in Sana et al. (2012), namely: the maxi- 
mum period at which significant binary interaction takes place is 
1500 d. Binaries with periods less than 6 days will initiate mass- 
transfer while the primary is still on the main sequence (Case A 
mass transfer). Among these systems, binaries with periods less 
than two days will merge. Between 2 and 6 days, binaries with 
mass-ratios less than 0.65 will also merge. Similarly a small frac- 
tion (5%) of case BC mass transfer (i.e., mass transfer occurring 
after the main sequence) will lead to coalescence. In the remain- 
ing interaction cases, the primary will be stripped from most/all 
of its envelope while the secondary will gain mass and angular 
momentum and will be spun up to critical rotational velocities. 

Under these assumptions and given an intrinsic binary frac- 
tion of 0.51 in the range 0.15 < log 10 P/d < 3.5, we estimate 
that 53% of all the stars born as O-type star belong to a binary 
system with a period less than 1500 d. The evolution of these 
stars is strongly affected by binary interaction: 18% of the O 
stars will merge with a companion, 27% will be stripped from 
their envelope and 8% are expected to be spun up. About one 
third of the binary interaction thus results in coalescence of the 
two companions. 

Compared to the frations of stars in the various binary evo- 
lution channels derived in the Milky Way (Sana et al. 2012), we 
obtain from the VFTS population an equivalent merger rate but 
lower rates of stripping and spin-up. In the above, we have im- 
plicitly assumed that today's distributions are representative of 
the distributions at birth. In an environment such as 30 Dor, this 
assumption may not hold. The presence of different populations 
in the field, some already quite old, and the fact that at least 5% 
of the O star population is running away, either because of dy- 
namical interaction or supernova kick, suggest that a fraction of 
the current O star population has already undergone significant 
dynamical and/or evolutionary interaction. 

In particular, today's binary fraction might be dragged to- 
wards lower values by the presence of runaways and post- 
interaction objects. As mentioned already the latter are indeed 
expected to be predominantly single (de Mink et al., in prepara- 
tion). Runaways and post-interaction objects would need to rep- 
resent about 20% of the current 30 Dor population to reconcile 
the binary fraction in 30 Dor (51%) with that of the Galactic 
open cluster population (69% Sana et al. 2012). Establishing 
the star formation and dynamical history of the Tarantula region 
would constitute a major step in deciding whether such a high 
fraction is plausible. 
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It remains that binary interaction effects have a critical influ- 
ence on the evolutionary path of more than half the stars born 
as O-type stars. These effects need to be considered to better in- 
terpret massive star populations seen in integrated light and to 
accurately explore the frequency of progenitors of compact ob- 
jects, high-mass X-ray binaries, hydrogen-poor supernovae and 
long-duration gamma-ray bursts. 

7. Conclusions 

In this paper, we investigate the multiplicity properties of a sam- 
ple of 360 O-type stars observed by multi-epoch spectroscopy 
in the framework of the VLT-Flames Tarantula Survey. Forty- 
six per cent of our sample present significant RV variations and, 
for 35%, the detected variations have an amplitude larger than 
20 kms -1 and are very likely spectroscopic binaries. 

The observed binary fraction present a high degree of homo- 
geneity across the field of view and among the various spectral- 
types and luminosity classes. Three sub-groups however show a 
significantly lower binary fraction: 

i. the binary fraction is lower in the outer region (r > 7.8'; 
=sll7 pc) of the field of view. This outer population may be 
dominated by single stars ejected from the core of the region; 

ii. the binary fraction is lower for the fainter stars in the sample 
(X s > 15.5), which results from the sample definition and 
traces the brightness separation between the O and B stars ; 

iii. the 09.7 I-II stars are predominantly single and belong to 
NGC 2060, the older O-star cluster in the field. Though spec- 
ulative at the moment, a large fraction of the 09.7 I-II stars 
may be post-interaction stars. 

We use a Monte-Carlo method to correct for the observa- 
tional biases and to constrain the intrinsic multiplicity proper- 
ties of the O-type star population. We obtain an intrinsic binary 
fraction of /bin = 0.5 1 ± 0.04. The most likely period distribu- 
tion, /(log 10 P/d) ~ (log 10 P/d)- 045 in the interval log 10 f/d e 
[0.15,3.5], favours shorter period systems compared to a flat 
distribution in log 10 P/d. Similarly the mass-ratio distribution, 
f(q) ~ q~ 10 in the interval q e [0.1, 1.0], favours lower-mass 
companions but we note that our method only provides weak 
constraints on the q distribution. 

The period distribution obtained here is strikingly similar 
to the one obtained for Galactic O-type binaries, possibly hint- 
ing a universal property among massive binaries. The intrinsic 
fraction of binaries with periods less that 10 3 5 d (0.51 ± 0.04) 
seems lower than obtained for the Galactic O-type binaries 
(0.69 ±0.09) albeit the two results agree within 2cr. Alternatively 
it could reflect the fact that the binary properties in the Tarantula 
region have already been significantly affected by dynamical 
and/or stellar evolution that would predominantly lead to either a 
merger event or disrupt the binary. Quantitatively understanding 
if and how binary evolution and dynamical interaction have af- 
fected the multiplicity properties in the Tarantula region would 
help to decide whether the observed differences are nature or 
nurture. A key ingredient in this process will be a better under- 
standing of the star formation history of the region. 

Our results emphasize again that multiplicity is one of the 
main characteristics of massive stars and that it needs to be taken 
into account to properly predict the evolution of entire popula- 
tions of massive stars. Dynamical interactions and evolutionary 
effects that could have possibly affected our sample would be ex- 
pected to, predominantly, decrease the observed number of bina- 
ries; our derived binary fraction can therefore be seen as a lower 



limit to the binary fraction at birth. This suggests that the most 
frequent final product of high-mass star formation is an O + OB 
binary and that massive star formation theories not only have to 
explain the formation of high mass stars but also have to repro- 
duce their multiplicity properties. 
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Appendix A: Normalisation 

In this section, we describe the semi-automatic procedure used 
to normalise the spectra discussed in this paper. The adopted al- 
gorithm uses all continuum points after automatically rejecting 
stellar and interstellar lines (both in emission and absorption) 
and cosmetic defects (bad columns, cosmic rays). 

The global trend in the spectrum is first removed using a one- 
degree polynomial fitted to the continuum, resulting in a spec- 
trum scaled around unity. A polynomial of degree 4 to 8 is then 
fitted to the scaled continuum. The error bars on the flux in each 
pixel are in principle provided by the ESO pipeline. By a com- 
parison of the error spectra with the empirical S /N in continuum 
regions, we find that the pipeline actually overestimates the er- 
ror spectra by a factor 1.2, which we correct for before using 
the pipeline-provided errors to weight the fit. The normalisation 
procedure proceeds as follows: 

i. the spectrum is smoothed using a 21 -pixel wide median fil- 
ter; 

ii. the spectral lines are removed by comparing the results of 
median and max/min filters applied on the initial spectrum, 
accounting for the spectrum S IN; 

iii. a first guess-solution of the continuum position is obtained 
using the smoothed spectrum (ii); 

iv. a second-guess solution is obtained using the smoothed spec- 
trum (ii) after cr-clipping around (iii); 

v. the final solution is obtained using the original spectrum (i.e., 
no smoothing) after cr-clipping around (iv). Manual exclu- 
sion of specific regions is allow at this stage to incorporate 
a priori knowledge on unsuitable continuum regions (e.g., 
broad diffuse interstellar bands, broad emission regions); 

vi. the spectrum (v) is scaled by the median value computed 
over its continuum regions. This allows to correct for slight 
global under/over estimates of the average continuum level 
which occurs when there are many absorption/emission lines 
as the line wings may induce a small but systematic bias. 
This empirical correction is of the order of a fraction of a per 
cent; 

vii. the normalisation function is applied to the error spectrum. 

Results are inspected. The degree of the polynomial and/or 
the exclusion of unsuitable continuum regions are modified 
until a satisfying fit is reached. We estimate that the continuum 
is, in general, constrained to better than 1% over the whole 
wavelength range, save for the edges of the spectra. 



Appendix B: Radial velocity measurements 

B.1. Methodology 

The choice of the RV measurement method is guided by the need 
to provide consistent, unbiased and homogeneous absolute RVs 
with representative error bars over the full range of O spectral 
sub-types and given the specificity of the data set. In our case, 
only a few lines are suitable for RV measurements. The number 
of epochs, the number of lines and their quality change from 
one star to another as a function of the star's spectral type and 
brightness and the degree of nebular contamination. The S /N 
also varies significantly from one star to another, from one epoch 
to another and over different spectral regions. 

Three RV techniques were considered: line moments, cross- 
correlation and Gaussian fitting. Tests were performed both on 



synthetic data and on a small set of objects taken from our sur- 
vey. Synthetic spectra for a range of temperatures, gravities, pro- 
jected rotational velocities (v sin z) and S /N representative of the 
VFTS O star data were generated using fastwind (Puis et al. 
2005) and degraded to the resolution of our survey. We briefly 
discuss the pros and cons of the three approaches below. 

Line moments turn out to be very sensitive to residuals of 
the nebular correction and the accuracy of the method is not ex- 
plored further. Auto cross-correlation (Dunstall et al., in prepa- 
ration) provides accurate RVs and is robust against residuals of 
the nebular correction. However, for stars with less than four 
lines suitable for RV measurement, cross-correlation has intrin- 
sic difficulties in providing accurate error bars. Gaussian fitting 
provides slightly less accurate RVs than cross-correlation, but 
with more reliable error bars for stars with only few lines of 
sufficient quality. However, the performance of Gaussian fitting 
degrades strongly with the quality of the lines (i.e., low S /N or 
high vsin/). 

To improve the performance of the Gaussian fitting, we de- 
veloped a tool to allow us to fit all the available lines at all epochs 
(and both wavelength set-ups) simultaneously. We assume that 
the line profiles are constant and well described by Gaussians. 
Each considered line is required to have the same full-width at 
half-maximum (FWHMs) and amplitude (A) at all epochs. All 
the lines of a given epoch are further required to display the same 
RV shift. If L is the number of considered lines and N the number 
of epochs, the number of parameters in the fit is thus : N + 2xL. 
Typical values for N are 18 and 6, depending on whether con- 
secutive exposures are taken individually or summed up. Values 
for L vary between 2 and 5. This yields a maximum number of 
parameters of the order of 28 in the present approach, to be com- 
pared with the 3 xWxL = 270 parameters of the standard ap- 
proach that adjusts each line separately. In essence the proposed 
method allows the line profile parameters (FWHMs and ampli- 
tudes) to be bootstrapped by the best quality spectra, improving 
the RV measurements of the lowest S/N data. 

For simplicity and robustness in the SB2 cases (see be- 
low), we have chosen to fit the full line profile using Gaussians. 
Intrinsic He profiles are however not always well represented by 
Gaussian shapes. For slow rotators (v rot sin/ < 80 kms~'), the 
shape of the line is dominated by the instrumental profiles that 
is well approximated by a Gaussian. The profile of moderately 
rotating stars (v rot sin; < 300 kms -1 ) is also well matched by a 
Gaussian, but with a small difference in the line wings. For fast 
rotators, rotationally broadened line profiles deviate significantly 
from the Gaussian shape. However, because we only attempt to 
measure the position of the centre of the lines and because the 
fitted lines are symmetric to within first-order, no significant bias 
is expected from this approximation. 

The fitting method is based on least-square estimates where 
the merit function accounts for the error spectrum to provide in- 
dividual error bars for each pixel. The optimisation relies on a 
modified Levensberq-Marquardt algorithm that requires a first 
estimate of the parameters. For constant RV stars and single- 
lined (SB1) binaries, the method is found to be robust against 
the initial guesses and we have taken v = 270 kms -1 , FWHM= 
3.0 A and A = 0.1 to initialize the fit. The allowed ranges for 
the various parameters are as follow: v e [-300, +700] kms -1 , 
A e [0.0,0.6] and FWHM e [0, 10] A. Only absorption lines 
were considered. For the doubled-lined spectroscopic binaries 
(SB2) in our sample, we follow the same approach with two 
Gaussian components per line profile. For some of the objects 
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Fig. B.l. Difference between the RVs measured from individual He lines and those measured from Hen /14541 as a function of the spectral sub- 
type. The panel at the bottom right corner compares RVs obtained using all available He lines (except He i+ii /14026 and He n /!4686) with the 
He ii /14541 RVs. The dashed line indicates the average of all comparisons weighted by their uncertainties. Only stars with constant RVs and data 
points with cr^v < 20 km s -1 have been included. 



H. Sana et al.: The VLT-FLAMES Tarantula Survey, Online Material p 3 




02 56 7 8 9 9.5 9.7 BO 
Stars sorted by sub-spectral type 



> 
I 



> 




02 56 7 8 9 9.5 9.7 BO 
Stars sorted by sub-spectral type 



Fig. B.2. Comparison of the RVs measured using simultaneously all suitable Hei and Hen lines (left and right panel, respectively) with the RVs 
measured using all suitable He lines. The Hei+n /14026 and He n /14686 lines have not been used in the measurements. 



the initial guesses needed to be iterated once or twice to improve 
the quality of the SB2 fit. 

B.2. Line selection 

In this section, we discuss the choice of spectral lines suitable 
for relative and absolute RV measurements. The spectra of O- 
type stars are dominated by Hi, Hei and Hen lines. Hydrogen 
lines are generally poorer RV indicators than (most of) the He 
lines because they are broader and more sensitive to wind ef- 
fects. The hydrogen lines of most of the VFTS O stars are also 
strongly affected by nebular emission. Metal lines (C, N, O, Si, 
Mg) are also present in O star spectra but are typically weaker 
than He lines. The presence of a given metal line is often lim- 
ited to a small range of spectral sub-types and does not allow for 
an homogeneous approach across the O star domain. Our anal- 
ysis is thus focused on the strong He i and He n lines from 4000 
to 5000 A, i.e. Hei AA43S7, 4471, 4713, 4922, Hen ,1,14200, 
4541, 4686 and the Hei+n blend at ,14026. In the rest of this 
section, we discuss the respective merits of each of these lines. 
Two aspects are important to consider: 

- the internal consistency, i.e. RV measurements from various 
lines of a given star should, coherently, reflect the systemic 
velocity of the star; 

- the homogeneity with respect to the full sample, i.e. the se- 
lected set of lines should provide RV measurements that can 
be compared to RVs from other stars even if only a sub- 
sample of the lines are available. This is important as the 
spectrum of early O stars is dominated by Hen lines and 
those of late O stars by He i lines. 

We first use simulated data to investigate the effect of tem- 
perature and wind strength on the measured RVs. Towards ear- 
lier spectral types, the He i+n blend at /14026 progressively shifts 
from pure Hei at ,14026.191 in late O stars to pure Hen at 
/14025.602 in early O stars. The RV measured from this blend 
can thus vary by over 30 kms -1 unless the effective rest wave- 
length of the blend is adapted as a function of the effective tem- 
perature of the star. Temperature also has an effect on the ac- 
curacy of the RV measurements. Hei lines provide more accu- 
rate RVs in mid to late O stars while He n lines are favoured for 
mid and early sub-types. The most accurate RVs are obtained for 
mid O-type stars because their spectra feature a larger number of 
good quality lines. 



Of all helium lines, stellar winds mostly affect He n ,14686, 
though not all spectral sub-types or luminosity class are equaaly 
affected by the filling in of the photospheric absorption line by 
wind emission. We find that RVs obtained from He n ,14686 are 
reliable for late and mid main-sequence sub-types. They are bi- 
ased towards negative values by up the 25 kms in early main- 
sequence stars. RVs from giants and supergiants are blue-shifted 
even for late-O stars and the line cannot be used as a probe for 
absolute RVs. 

We also compare RV measurements obtained individually 
from different lines in a representative sub-sample of the O stars, 
covering the complete O spectral sub-range. Only stars display- 
ing constant RVs are considered for this test. Because the line 
is observed both in the LR02 and LR03 setups and remains of 
a reasonable strength even in the latest O sub-types, we use 
the He ii ,14541 line as a reference for the comparison. Fig. B.l 
shows the difference between RVs obtained for a given line and 
those obtained from He n ,14541 plotted against the spectral sub- 
type. Temperature effects on He i+ii ,14026 are clearly illustrated. 
Wind effects turn out to have a limited impact on Hen ,14686 
measurements mostly because our sample is dominated by main- 
sequence stars and because many single supergiants present low- 
amplitude RV variations and are thus excluded from the compar- 
ison. 

RVs from the Hei ,1447 1 line show large discrepancies, 
which we attribute to the combination of the following factors: 
He i ,14471 is a triplet transition, with a forbidden component. It 
is quasi-metastable, similar to but less extreme than He i A/13889, 
5876, 10830 and is highly susceptible to wind effects. For late 
spectral sub-types it is blended with On in its blue wing and, 
among the measured He i lines, it suffers the most from the neb- 
ular contamination. 

For late spectral types, Hen ,14200 is blended with 
Nra ,14200.07 (i.e., about 17 km/s redwards from Hei ,14200), 
which can explain some of the deviations seen in the correspond- 
ing panel in Fig. B.l. For late O stars, the uncertainties associ- 
ated with He ii ,14200 RVs are large anyway and the line has a 
very limited weight in the final RV measurements obtained from 
simultaneously fitting all available lines (Fig. B.2). 

The best consistency is observed between the lines 
He ii ,1,14200, 4541, He i ,1,14387 and 4713. He i ,14922 displays 
good consistency over the range of O spectral sub-types but pro- 
vides RVs that are on average shifted by 3 kms -1 compared to 
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Table B.l. Rest wavelengths (A Q ) used to computed the RVs. 



Line 


<lo (A) 


Reference 


1 i_i_n 2A(Y)f> 
rlc I+II A^+KJZO 


<+UZO.U / Z 


CJ T 77 


He ii ^4200 


4199.832 


PvH/NIST 


He 1,14387 


4387.929 


PvH 


He 1,14471 


4471.480 


NIST 


He ii ,14541 


4541.591 


PvH/NIST 


He ii ,14686 


4685.710 


PvH 


Hei ,14713 


4713.146 


NIST 


He 1/14922 


4921.931 


PvH/NIST 



References: CLL77: Conti et al. (1977); NIST atomic spectra data 
base: http://www.nist.gov/pml/data/asd.cfm; PvH : Peter van 
Hoof atomic line list (http://www.pa.uky.edu/~peter/atomic/). 
Notes: a. Because of the strong blend between Hei /14026.1914 and 
He ii ,14025.602, the effective rest wavelength is function of the 
temperature, thus the spectral type. We use here an observed value 
which is representative of mid- to late- O-type stars. 

the other lines. We mitigate this by modifying the effective rest 
wavelength of the transition. 

As a conclusion, up to five lines can be used for absolute 
RV measurements (He n /14200, 4541, He i /U4387, 4713, 4922) 
depending on the spectral types, S /N and nebular contamination. 
He i+ii /14026, and, for relatively weak wind stars, He n ,14686 
can be used for relative measurements. 

Appendix C: A Monte-Carlo method to constrain 
the intrinsic orbital parameter distributions 

In this section, we provide the details of the Monte-Carlo method 
used to constrain the intrinsic orbital parameter distributions. 
The general philosophy resembles that of Kobulnicky & Fryer 
(2007), in the sense that it relies on a Monte-Carlo modeling of 
the orbital velocities. On the one hand, we follow Kobulnicky 
& Fryer in using KS tests to compare simulated and observed 
distributions and the diagnostic distributions based on ARV and 
X 1 (See Sect. C.3) are similar to, respectively, the V h and V rms 
used by Kobulnicky & Fryer. On the other hand, our method 
fundamentally differs in various aspects. First we use the specific 
measurement errors provided by the RV analysis throughout the 
complete Monte-Carlo analysis. Second, we do not apply exter- 
nal constraints such as fixing the period distribution or using the 
fraction of Type Ib/c supernovae corresponding to each model. 
We thus attempt to simultaneously constrain the intrinsic period 
and mass-ratio distributions and the intrinsic binary fraction. 

The core of the method is identical to that described in Sana 
et al. (2012) but differs in the adopted merit function. Sana et al. 
compared the observed orbital parameter distributions while we 
will use distributions constructed from the RV database. Indeed 
most of the detected binaries in 30 Dor have too few RV mea- 
surements to compute a meaningful orbital solution, precluding 
the direct approach used by Sana et al. (2012). 

C.1. Method overview 

As discussed in the main text, we use power laws to describe 
the intrinsic distributions of orbital parameters: /(log 10 P/d) ~ 
(log^P/dy 1 , f(q) ~ q" and /(e) ~ q n , the exponents of which 
are left as free parameters. For each combination of n, k and /bin 
in our grid, we draw populations of /V = 360 primary stars using 
a Kroupa mass function (Kroupa 2001) between 15 and 80 M Q , 
covering thus the expected mass range of O-type stars. Each star 



is assigned an observational sequence (observing epochs and RV 
accuracy at each epoch) from our VFTS sample. A fraction of 
these stars are randomly assigned to be binaries following a bi- 
nomial with a success probability /,;„. The orbital parameters of 
the binaries are randomly drawn as follows : 

- the period and mass-ratio are taken from the power-law dis- 
tributions given 7r and k; 

- the eccentricity is taken from a power-law distribution with 
an index given by r\ — -0.5 (see discussion in Sect. 4.1); 

- the inclination and periastron argument are taken assuming 
random orientation of the orbit in space; 

- the time of periastron passage is assumed to be uncorrelated 
with the observational sampling and is thus taken from a uni- 
form distribution. 

For each binary, the orbital velocities associated with the pri- 
mary stars are computed taking into account the time sequence 
of our observations. Measurement uncertainties are randomly 
added using a Gaussian distribution with a standard deviation 
given by the associated measurement error. For simulated single 
stars, only the contribution of the error is taken into account. 

The binary detection criteria of Eq. 4 are then applied and 
the simulated detected binaries are flagged. Simulated observed 
distributions (see Sect. C.3) are built based on the simulated RV 
measurements of the detected systems and the simulated mea- 
sured binary fraction /J™" 1 is computed. 

For a given combination of n, k and /bin the process is re- 
peated 100 times to build the parent statistics. The simulated 
parent distributions are then compared to the observed ones by 
means of a one-sided KS test. We also compute the binomial 
probability B(Nbi n ,N, Z^™" 1 ) to detect Nu n binaries among a pop- 
ulation of N stars given the success probability predicted by the 
simulations, i.e. /Sf" 1 . Adopting Nu n to be the actual number 
of detected O-type binaries in VFTS allows us to estimate the 
ability of the assumed intrinsic parameters to reproduce the ob- 
served binary fraction while taking into account the finite size of 
the sample. 

C.2. Diagnostic distributions 

In this section, we specifically discuss the empirical distributions 
used as diagnostics. The choice of the merit function will be dis- 
cussed in the next section. We investigate four diagnostic quan- 
tities: 

i. the distribution of, for each detected binary, the amplitude of 
the RV signal (ARV); 

ii. the distribution of, for each detected binary, the^f 2 of a con- 
stant RV model; 

iii. the distribution of, for each detected binary, the minimum 
time scale for significant RV variations to be seen (AHJD); 

iv. the detected binary fraction. 

We thus perform a parameter study to investigate the effect 
of different exponents on the simulated distributions and on the 
overall binary detection rate. Figure C.l shows that the exponent 
n of the period distribution has the largest impact on the simu- 
lated distributions. The exponent k of the mass-ratio distribution 
is of secondary importance and the exponent rj has almost no 
impact on the simulated distribution, albeit affecting the simu- 
lated detection rate by a couple of per cent. The distribution of 
the minimum time scale for the observed RV variations is al- 
most exclusively sensitive to the period distribution but mostly 
independent of the mass-ratio. 
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Fig. C.l. Comparison between the observed (crosses) and simulated (lines) cumulative distributions of the peak-to-peak; RV amplitudes (top row), 
variability time scales (middle row) and x 2 (bottom row). Lefthand panels vary the exponent n of the period distribution. Panels in the central 
column vary k, the exponent of the mass-ratio distribution and the righthand panels vary 77, the exponent of the eccentricity distribution, n, k and 
77 values are indicated in the upper-left legend of each panel. The bottom-right legends indicate the overall VFTS detection probability for the 
adopted parent distributions. 



C.3. Suitable merit functions 

At each point of the grid defined by the the intrinsic values n, 
k and /bi n , the simulated distributions (i) to (iii) are individually 
compared to the observed ones using KS-tests while the sim- 
ulated and observed binary fraction are compared using a bi- 
nomial statistics. Because we aim at matching all the diagnos- 
tic distributions simultaneously, we also combine the individual 
probabilities in a global merit function. 



We run several test to assess different merit functions (i.e., 
different ways of combining the individual probabilities) as well 
as to validate our method. We generate synthetic data from 
known input distributions, i.e. known n, k and /b; n , and we apply 
our code using different merit functions. Input and output distri- 
butions are then compared to check the suitability of the various 
merit functions as well as the general accuracy of the method. 
The process has been repeated using 50 synthetic data sets in 
each case. The synthetic data have been generated such that they 
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Table C.l. Overview of the test results using different merit functions (Col. 2). The median, 0.16 and 0.74 percentiles of the retrieved parameters 
in a set of 50 test runs are given in Cols. 4 to 6. The input values of these parameters are indicated as the second line of the table. The explicit 
mathematical descriptions are given for the individual approaches (a)-(c) and (l)-(m). Implicit abbreviations are used otherwise to keep the notation 
compact. 
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share the same observational properties (sampling, measurement 
accuracy) as the VFTS Medusa data. Table C.l provides an 
overview of the results of the tests. Because of the computa- 
tional load of such tests, we use a slightly coarser grid with steps 
of 0.02 in /bin- Explored ranges in /bin, x and k and steps in 
the power law indexes n and k are identical to those quoted in 



Table 5. In addition to appraising the merit function, we also test 
the need to apply a minimal threshold C for the detected signifi- 
cant RV signal. 

The first striking result is that most merit functions can re- 
cover the input binary fraction within a few percents albeit with 
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various precision. Some of the best results are obtained when 
including the B(N b s , N, / S i m ) term in the merit function. 

The index n of the period distribution is poorly constrained 
whenever one uses no detection threshold (i.e., C = kms 1 ). 
This indicates that the detection threshold is not only useful to 
reject false detections due to intrinsic profile variability but is 
also a critical ingredient of the method. Interestingly, n is reason- 
ably well constrained by the AHJD distribution alone, i.e. merit 
function (1), but a better overall result is obtained when used in 
conjunction with other diagnostic quantities. 

The ARV and^- 2 distributions are useful to refine the k value 
but the uncertainties remain large. In essence, the ARV and x 2 
distributions carry similar information. Because merit functions 
using ARV tend to behave slightly better than those using x 2 , 
and because the x 2 values are by definition much more sensitive 
to the exact knowledge of the measurement errors, we select our 
final merit function as given by the product of the Pks probabil- 
ities computed from the AHJD and ARV distributions and of the 
Binomial probability (Eq. 5), hence merit function (u) in Table 
C.l. 

Overall, the retrieved values of n and k tend to be smaller 
than their input values but the correct indexes concur with the 
lcr confidence intervals quoted in Table C.l. These confidence 
intervals can further be used to estimate the formal uncertainty of 
the method. Given that the best representation of the data are ob- 
tained with / b i n = 0.51, n = -0.45 and k = -1.00 (see Sect. 4.2), 
the bottom part of Table C.l applies and we adopt o-f bin = 0.04, 
cr n = 0.3 and cr K = 0.4. 



Appendix D: Supplementary figures 
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Fig. D.l. A^ s -band magnitudes as a function of the spectral sub-types for luminosity classes V-IV (left panel), III (middle panel) and II-I (right 
panel). Different symbols identify the single stars, the low-amplitude RV variables and the binaries. Symbols for single and binaries of a given 
spectral-types have been shifted, compared to one another, by a small amount amount along the x-axis to improve the clarity of the plots. 
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Fig. D.2. Same as Fig. 1 1 for the luminosity classes V-IV (left panels), III (middle panels) and II-I (right panels). 
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Fig.D.3. Same as Fig. 12 for the luminosity classes V-IV (left panels), III (middle panels) and II-I (right panels). 



